Vowel Acoustics as Predictors of Speech Intelligibility in Dysarthria

This is the code for the statistical analysis for “Vowel Acoustics as Predictors of Speech Intelligibility in Dysarthria.”

Loading Packages

library(rio)

The following rio suggested packages are not installed: ‘csvy’, ‘feather’, ‘fst’, ‘hexView’, ‘readODS’, ‘rmatio’
Use 'install_formats()' to install them

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ─────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.5     ✓ purrr   0.3.4
✓ tibble  3.1.6     ✓ dplyr   1.0.7
✓ tidyr   1.1.3     ✓ stringr 1.4.0
✓ readr   1.4.0     ✓ forcats 0.5.1
── Conflicts ────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(irr) # install.packages('irr')

Loading required package: lpSolve

library(performance)

Attaching package: ‘performance’

The following object is masked from ‘package:irr’:

    icc

library(car)

Loading required package: carData

Attaching package: ‘car’

The following object is masked from ‘package:dplyr’:

    recode

The following object is masked from ‘package:purrr’:

    some

library(ggpubr)
library("Hmisc") # install.packages('Hmisc')

Loading required package: lattice
Loading required package: survival
Loading required package: Formula

Attaching package: ‘Hmisc’

The following objects are masked from ‘package:dplyr’:

    src, summarize

The following objects are masked from ‘package:base’:

    format.pval, units

library(ggridges)
library(furniture) # install.packages('furniture')
library(gt)

Attaching package: ‘gt’

The following object is masked from ‘package:Hmisc’:

    html

library(patchwork)
library(ks)
library(emuR) # install.packages('emuR')

Attaching package: ‘emuR’

The following object is masked from ‘package:Hmisc’:

    label

The following object is masked from ‘package:car’:

    ellipse

The following object is masked from ‘package:base’:

    norm

Upload Datasets

Reliability <- rio::import("Prepped Data/Reliability Data.csv")
AcousticData <- rio::import("Prepped Data/AcousticMeasures.csv") %>%
  dplyr::mutate(intDiff = VAS - transAcc)

AcousticData <- AcousticData %>%
  dplyr::filter(!grepl("_rel", Speaker)) %>%
  dplyr::select(c(Speaker, Sex, Etiology, vowel_ED_b, VSA_b, autoVSA,
                  Hull_b,Hull_bVSD_25, Hull_bVSD_50, Hull_bVSD_75,
                  VAS, transAcc)) %>%
  dplyr::mutate(Etiology = as.factor(Etiology),
                Sex = as.factor(Sex),
                Speaker = as.factor(Speaker))

Listeners <- rio::import("Prepped Data/Listener_Demographics.csv") %>%
  dplyr::select(!c(StartDate:proloficID, Q2.4_6_TEXT, Q3.2_8_TEXT, AudioCheck:EP3))

Listeners$race[Listeners$Q3.3_7_TEXT == "Native American/ African amercing"] <- "Biracial or Multiracial"

Inter-rater Reliability

Two raters (the first two authors) completed vowel segmentation for the speakers. To calculate inter-rater reliability, 20% of the speakers were segmented again by the other rater. Two-way intraclass coefficients were computed for the extracted F1 and F2 from the temporal midpoint of the vowel segments. Since only one set of ratings will be used in the data analysis, we focus on the single ICC results and interpretation. However, we also report the average ICC values to be comprehensive.

## Creating new data frames to calculate ICC for extracted F1 and F2 values

F1_Rel <- Reliability %>%
  dplyr::select(c(F1, F1_rel))

F2_Rel <- Reliability %>%
  dplyr::select(c(F2, F2_rel))
  
## Single ICC for F1
Single_F1 <- irr::icc(F1_Rel, model = "twoway", type = "agreement", unit = "single")

## Average ICC for F1
Average_F1 <- irr::icc(F1_Rel, model = "twoway", type = "agreement", unit = "average")

## Single ICC for F2
Single_F2 <- irr::icc(F2_Rel, model = "twoway", type = "agreement", unit = "single")

## Average ICC for F2
Average_F2 <- irr::icc(F2_Rel, model = "twoway", type = "agreement", unit = "average")

## Inter-rater reliability results and interpretation

print(paste("Single ICC for F1 is ", round(Single_F1$value, digits = 3), ". ", 
            "The 95% CI is [", round(Single_F1$lbound, digits = 3), " - ", round(Single_F1$ubound, digits = 3), "].", sep = ""))

[1] "Single ICC for F1 is 0.866. The 95% CI is [0.837 - 0.89]."

print(paste("Single ICC for F2 is ", round(Single_F2$value, digits = 3), ". ", 
            "The 95% CI is [", round(Single_F2$lbound, digits = 3), " - ", round(Single_F2$ubound, digits = 3), "].", sep = ""))

[1] "Single ICC for F2 is 0.931. The 95% CI is [0.916 - 0.944]."

print(paste("Average ICC for F1 is ", round(Average_F1$value, digits = 3), ". ", 
            "The 95% CI is [", round(Average_F1$lbound, digits = 3), " - ", round(Average_F1$ubound, digits = 3), "].", sep = ""))

[1] "Average ICC for F1 is 0.928. The 95% CI is [0.911 - 0.942]."

print(paste("Average ICC for F2 is ", round(Average_F2$value, digits = 3), ". ", 
            "The 95% CI is [", round(Average_F2$lbound, digits = 3), " - ", round(Average_F2$ubound, digits = 3), "].", sep = ""))

[1] "Average ICC for F2 is 0.964. The 95% CI is [0.956 - 0.971]."

print("Thus, interrater reliability for the extracted F1 and F2 values from the vowel segments was good to excellent.")

[1] "Thus, interrater reliability for the extracted F1 and F2 values from the vowel segments was good to excellent."

## Removing extra data frames from environment

rm(F1_Rel, F2_Rel, Reliability, Single_F1, Single_F2, Average_F1, Average_F2)

Descriptive Statistics

Correlations

CorrMatrix <- AcousticData %>%
  dplyr::select(VSA_b, vowel_ED_b, autoVSA, Hull_b, Hull_bVSD_25, Hull_bVSD_75, VAS, transAcc) %>%
  as.matrix() %>%
  Hmisc::rcorr()

CorrMatrixP <- CorrMatrix$P < .05

CorrMatrix <- CorrMatrix$r

stats::cor.test(AcousticData$VSA_b, AcousticData$vowel_ED_b, method = "pearson")

    Pearson's product-moment correlation

data:  AcousticData$VSA_b and AcousticData$vowel_ED_b
t = 6.5285, df = 38, p-value = 1.076e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5372669 0.8468014
sample estimates:
      cor 
0.7270881 

stats::cor.test(AcousticData$Hull_b, AcousticData$Hull_bVSD_25, method = "pearson")

    Pearson's product-moment correlation

data:  AcousticData$Hull_b and AcousticData$Hull_bVSD_25
t = 9.4906, df = 38, p-value = 1.43e-11
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.7135136 0.9119080
sample estimates:
     cor 
0.838625 

stats::cor.test(AcousticData$Hull_b, AcousticData$autoVSA, method = "pearson")

    Pearson's product-moment correlation

data:  AcousticData$Hull_b and AcousticData$autoVSA
t = 11.103, df = 38, p-value = 1.723e-13
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.7734205 0.9319757
sample estimates:
      cor 
0.8742893 

write.csv(CorrMatrix, file = "Tables/Correlation Matrix.csv")
rm(CorrMatrix)

Research Q1: Modeling Intelligibility

Orthographic Transcriptions

Model 1

# Specifying Model 1

OT_Model1 <- lm(transAcc ~ Hull_bVSD_25, data = AcousticData)

## Model 1 Assumptions 

performance::check_model(OT_Model1)

## Model 1 Summary

summary(OT_Model1)

Call:
lm(formula = transAcc ~ Hull_bVSD_25, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-47.896 -14.090   5.996  17.470  36.105 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   48.2973     9.7372   4.960  1.5e-05 ***
Hull_bVSD_25   0.6417     0.5663   1.133    0.264    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 23.83 on 38 degrees of freedom
Multiple R-squared:  0.03268,   Adjusted R-squared:  0.007224 
F-statistic: 1.284 on 1 and 38 DF,  p-value: 0.2643

Model 2

## Specifying Model 2

OT_Model2 <- lm(transAcc ~ Hull_bVSD_25 + Hull_bVSD_75, data = AcousticData)

## Model 2 Assumption Check

performance::check_model(OT_Model2)

## Model 2 Summary

summary(OT_Model2)

Call:
lm(formula = transAcc ~ Hull_bVSD_25 + Hull_bVSD_75, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-48.854 -13.962   5.567  16.758  36.257 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   47.3374    10.3316   4.582 5.09e-05 ***
Hull_bVSD_25   0.8045     0.7781   1.034    0.308    
Hull_bVSD_75  -0.7188     2.3225  -0.309    0.759    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 24.11 on 37 degrees of freedom
Multiple R-squared:  0.03518,   Adjusted R-squared:  -0.01698 
F-statistic: 0.6745 on 2 and 37 DF,  p-value: 0.5156

## Model 1 and Model 2 Comparison

anova(OT_Model1, OT_Model2)

Analysis of Variance Table

Model 1: transAcc ~ Hull_bVSD_25
Model 2: transAcc ~ Hull_bVSD_25 + Hull_bVSD_75
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1     38 21570                           
2     37 21514  1    55.696 0.0958 0.7587

Model 3a

## Specifying Model 3

OT_Model3a <- lm(transAcc ~ Hull_bVSD_25 + Hull_bVSD_75 + Hull_b, data = AcousticData)

## Model 3 Assumption Check

performance::check_model(OT_Model3a)
performance::check_collinearity(OT_Model3a)


## Model 3 Summary

summary(OT_Model3a)

## Model 2 and Model 3 Comparison

anova(OT_Model2, OT_Model3a)

Model 3b

## Specifying Model 3

OT_Model3b <- lm(transAcc ~ Hull_bVSD_75 + Hull_b, data = AcousticData)

## Model 3 Assumption Check

performance::check_model(OT_Model3b)

performance::check_collinearity(OT_Model3b)

# Check for Multicollinearity

Low Correlation

         Term  VIF Increased SE Tolerance
 Hull_bVSD_75 1.26         1.12      0.79
       Hull_b 1.26         1.12      0.79

## Model 3 Summary

summary(OT_Model3b)

Call:
lm(formula = transAcc ~ Hull_bVSD_75 + Hull_b, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-54.371 -12.860   5.038  17.725  31.609 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept)   33.2337    13.9973   2.374   0.0229 *
Hull_bVSD_75  -0.6193     1.8702  -0.331   0.7424  
Hull_b         0.8605     0.4809   1.789   0.0818 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 23.47 on 37 degrees of freedom
Multiple R-squared:  0.08635,   Adjusted R-squared:  0.03697 
F-statistic: 1.749 on 2 and 37 DF,  p-value: 0.1881

Model 4a

## Specifying Model 4

OT_Model4a <- lm(transAcc ~ Hull_bVSD_75 + Hull_b + autoVSA, data = AcousticData)

## Model 4 Assumption Check

performance::check_model(OT_Model4a)

performance::check_collinearity(OT_Model4a)

# Check for Multicollinearity

Low Correlation

         Term  VIF Increased SE Tolerance
 Hull_bVSD_75 1.31         1.15      0.76
      autoVSA 4.41         2.10      0.23

Moderate Correlation

   Term  VIF Increased SE Tolerance
 Hull_b 5.02         2.24      0.20

## Model 4 Summary

summary(OT_Model4a)

Call:
lm(formula = transAcc ~ Hull_bVSD_75 + Hull_b + autoVSA, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-46.729  -9.736   2.928  13.425  32.586 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)   33.0695    13.0702   2.530  0.01592 * 
Hull_bVSD_75  -1.5011     1.7806  -0.843  0.40478   
Hull_b         2.8262     0.8956   3.156  0.00323 **
autoVSA       -5.6340     2.2208  -2.537  0.01566 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21.91 on 36 degrees of freedom
Multiple R-squared:  0.2249,    Adjusted R-squared:  0.1603 
F-statistic: 3.482 on 3 and 36 DF,  p-value: 0.02561

## Model 3 and Model 4 Comparison

anova(OT_Model3b, OT_Model4a)

Analysis of Variance Table

Model 1: transAcc ~ Hull_bVSD_75 + Hull_b
Model 2: transAcc ~ Hull_bVSD_75 + Hull_b + autoVSA
  Res.Df   RSS Df Sum of Sq      F  Pr(>F)  
1     37 20373                              
2     36 17283  1    3089.9 6.4361 0.01566 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Model 4b

## Specifying Model 4

OT_Model4b <- lm(transAcc ~ Hull_bVSD_75 + autoVSA, data = AcousticData)

## Model 4 Assumption Check

performance::check_model(OT_Model4b)

performance::check_collinearity(OT_Model4b)

# Check for Multicollinearity

Low Correlation

         Term  VIF Increased SE Tolerance
 Hull_bVSD_75 1.11         1.05      0.90
      autoVSA 1.11         1.05      0.90

## Model 4 Summary

summary(OT_Model4b)

Call:
lm(formula = transAcc ~ Hull_bVSD_75 + autoVSA, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-49.900 -15.387   6.626  17.304  32.102 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   52.3999    12.8682   4.072 0.000236 ***
Hull_bVSD_75   0.7069     1.8248   0.387 0.700710    
autoVSA        0.4296     1.2411   0.346 0.731176    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 24.42 on 37 degrees of freedom
Multiple R-squared:  0.0105,    Adjusted R-squared:  -0.04298 
F-statistic: 0.1964 on 2 and 37 DF,  p-value: 0.8225

Model 5

## Specifying Model 5

OT_Model5 <- lm(transAcc ~ Hull_bVSD_75 + autoVSA + VSA_b, data = AcousticData)

## Model 4 Assumption Check

performance::check_model(OT_Model5)

performance::check_collinearity(OT_Model5)

# Check for Multicollinearity

Low Correlation

         Term  VIF Increased SE Tolerance
 Hull_bVSD_75 1.15         1.07      0.87
      autoVSA 1.20         1.09      0.84
        VSA_b 1.17         1.08      0.85

## Model 4 Summary

summary(OT_Model5)

Call:
lm(formula = transAcc ~ Hull_bVSD_75 + autoVSA + VSA_b, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-46.674 -12.643   4.169  15.099  34.515 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)   38.1376    11.9031   3.204  0.00284 **
Hull_bVSD_75  -0.4416     1.6222  -0.272  0.78701   
autoVSA       -0.6501     1.1229  -0.579  0.56622   
VSA_b          6.3560     1.7810   3.569  0.00104 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21.28 on 36 degrees of freedom
Multiple R-squared:  0.2691,    Adjusted R-squared:  0.2082 
F-statistic: 4.418 on 3 and 36 DF,  p-value: 0.009589

## Model 3 and Model 4 Comparison

anova(OT_Model4b, OT_Model5)

Analysis of Variance Table

Model 1: transAcc ~ Hull_bVSD_75 + autoVSA
Model 2: transAcc ~ Hull_bVSD_75 + autoVSA + VSA_b
  Res.Df   RSS Df Sum of Sq      F   Pr(>F)   
1     37 22064                                
2     36 16298  1    5766.1 12.736 0.001038 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Model 6

## Specifying Model 6

OT_Model6 <- lm(transAcc ~ Hull_bVSD_75 + autoVSA + VSA_b + vowel_ED_b, data = AcousticData)

## Model 6 Assumption Check

performance::check_model(OT_Model6)

performance::check_collinearity(OT_Model6)

# Check for Multicollinearity

Low Correlation

         Term  VIF Increased SE Tolerance
 Hull_bVSD_75 1.20         1.10      0.83
      autoVSA 1.31         1.15      0.76
        VSA_b 2.30         1.52      0.43
   vowel_ED_b 2.36         1.54      0.42

## Model 6 Summary

summary(OT_Model6)

Call:
lm(formula = transAcc ~ Hull_bVSD_75 + autoVSA + VSA_b + vowel_ED_b, 
    data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-45.023 -13.342   2.208  17.093  32.847 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept)   26.9818    20.5367   1.314   0.1974  
Hull_bVSD_75  -0.2272     1.6659  -0.136   0.8923  
autoVSA       -0.8851     1.1848  -0.747   0.4600  
VSA_b          5.1767     2.5152   2.058   0.0471 *
vowel_ED_b     8.9714    13.4053   0.669   0.5077  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21.44 on 35 degrees of freedom
Multiple R-squared:  0.2783,    Adjusted R-squared:  0.1958 
F-statistic: 3.375 on 4 and 35 DF,  p-value: 0.01948

## Model 5 and Model 6 Comparison

anova(OT_Model5, OT_Model6)

Analysis of Variance Table

Model 1: transAcc ~ Hull_bVSD_75 + autoVSA + VSA_b
Model 2: transAcc ~ Hull_bVSD_75 + autoVSA + VSA_b + vowel_ED_b
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1     36 16298                           
2     35 16092  1    205.93 0.4479 0.5077

Final Model

## Specifying Final Model

OT_Model_final <- lm(transAcc ~ VSA_b, data = AcousticData)

## Final Model Assumption Check

performance::check_model(OT_Model_final)

## Final Model Summary

summary(OT_Model_final)

Call:
lm(formula = transAcc ~ VSA_b, data = AcousticData)

Residuals:
   Min     1Q Median     3Q    Max 
-46.72 -12.69   2.97  14.37  35.39 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   32.508      7.857   4.138 0.000187 ***
VSA_b          5.872      1.613   3.641 0.000807 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 20.86 on 38 degrees of freedom
Multiple R-squared:  0.2586,    Adjusted R-squared:  0.2391 
F-statistic: 13.25 on 1 and 38 DF,  p-value: 0.0008068

confint(OT_Model_final)

                2.5 %    97.5 %
(Intercept) 16.602765 48.413532
VSA_b        2.606927  9.137097

VAS Models

Model 1

# Specifying Model 1

VAS_Model1 <- lm(VAS ~ Hull_bVSD_25, data = AcousticData)

## Model 1 Assumptions 

performance::check_model(VAS_Model1)

## Model 1 Summary

summary(VAS_Model1)

Model 2

## Specifying Model 2

VAS_Model2 <- lm(VAS ~ Hull_bVSD_25 + Hull_bVSD_75, data = AcousticData)

## Model 2 Assumption Check

performance::check_model(VAS_Model2)

## Model 2 Summary

summary(VAS_Model2)

Call:
lm(formula = VAS ~ Hull_bVSD_25 + Hull_bVSD_75, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-47.850 -16.576   8.382  19.448  37.237 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   42.3384    11.4211   3.707 0.000684 ***
Hull_bVSD_25   0.6376     0.8602   0.741 0.463195    
Hull_bVSD_75  -0.2204     2.5674  -0.086 0.932036    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 26.66 on 37 degrees of freedom
Multiple R-squared:  0.02291,   Adjusted R-squared:  -0.0299 
F-statistic: 0.4338 on 2 and 37 DF,  p-value: 0.6513

## Model 1 and Model 2 Comparison

anova(VAS_Model1, VAS_Model2)

Analysis of Variance Table

Model 1: VAS ~ Hull_bVSD_25
Model 2: VAS ~ Hull_bVSD_25 + Hull_bVSD_75
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1     38 26296                           
2     37 26291  1     5.239 0.0074  0.932

Model 3a

## Specifying Model 3

VAS_Model3a <- lm(VAS ~ Hull_bVSD_25 + Hull_bVSD_75 + Hull_b, data = AcousticData)

## Model 3 Assumption Check

performance::check_model(VAS_Model3a)

performance::check_collinearity(VAS_Model3a)

# Check for Multicollinearity

Low Correlation

         Term  VIF Increased SE Tolerance
 Hull_bVSD_75 2.00         1.41      0.50
       Hull_b 3.65         1.91      0.27

Moderate Correlation

         Term  VIF Increased SE Tolerance
 Hull_bVSD_25 5.33         2.31      0.19

## Model 3 Summary

summary(VAS_Model3a)

Call:
lm(formula = VAS ~ Hull_bVSD_25 + Hull_bVSD_75 + Hull_b, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-56.444 -18.989   6.121  18.036  32.281 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)   25.0400    16.0599   1.559    0.128
Hull_bVSD_25  -1.1164     1.4389  -0.776    0.443
Hull_bVSD_75   0.8805     2.6280   0.335    0.740
Hull_b         1.3770     0.9139   1.507    0.141

Residual standard error: 26.21 on 36 degrees of freedom
Multiple R-squared:  0.08087,   Adjusted R-squared:  0.00428 
F-statistic: 1.056 on 3 and 36 DF,  p-value: 0.3799

## Model 2 and Model 3 Comparison

anova(VAS_Model2, VAS_Model3a)

Analysis of Variance Table

Model 1: VAS ~ Hull_bVSD_25 + Hull_bVSD_75
Model 2: VAS ~ Hull_bVSD_25 + Hull_bVSD_75 + Hull_b
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1     37 26291                           
2     36 24731  1    1559.6 2.2702 0.1406

Model 3b

## Specifying Model 3b

VAS_Model3b <- lm(VAS ~  Hull_bVSD_75 + Hull_b, data = AcousticData)

## Model 3 Assumption Check

performance::check_model(VAS_Model3b)

performance::check_collinearity(VAS_Model3b)

# Check for Multicollinearity

Low Correlation

         Term  VIF Increased SE Tolerance
 Hull_bVSD_75 1.26         1.12      0.79
       Hull_b 1.26         1.12      0.79

## Model 3 Summary

summary(VAS_Model3b)

Call:
lm(formula = VAS ~ Hull_bVSD_75 + Hull_b, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-53.356 -19.394   8.036  21.298  33.412 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept)   27.8884    15.5503   1.793   0.0811 .
Hull_bVSD_75  -0.3567     2.0777  -0.172   0.8646  
Hull_b         0.8034     0.5343   1.504   0.1412  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 26.07 on 37 degrees of freedom
Multiple R-squared:  0.0655,    Adjusted R-squared:  0.01499 
F-statistic: 1.297 on 2 and 37 DF,  p-value: 0.2855

Model 4a

## Specifying Model 4

VAS_Model4a <- lm(VAS ~ Hull_bVSD_75 + autoVSA + Hull_b, data = AcousticData)

## Model 4 Assumption Check

performance::check_model(VAS_Model4a)

performance::check_collinearity(VAS_Model4a)

# Check for Multicollinearity

Low Correlation

         Term  VIF Increased SE Tolerance
 Hull_bVSD_75 1.31         1.15      0.76
      autoVSA 4.41         2.10      0.23

Moderate Correlation

   Term  VIF Increased SE Tolerance
 Hull_b 5.02         2.24      0.20

## Model 4 Summary

summary(VAS_Model4a)

Call:
lm(formula = VAS ~ Hull_bVSD_75 + autoVSA + Hull_b, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-51.912 -20.901   3.669  16.070  40.028 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)   27.7029    14.4764   1.914  0.06364 . 
Hull_bVSD_75  -1.3527     1.9721  -0.686  0.49717   
autoVSA       -6.3639     2.4597  -2.587  0.01386 * 
Hull_b         3.0238     0.9919   3.048  0.00429 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 24.27 on 36 degrees of freedom
Multiple R-squared:  0.212, Adjusted R-squared:  0.1464 
F-statistic: 3.229 on 3 and 36 DF,  p-value: 0.03365

## Model 3 and Model 4 Comparison

anova(VAS_Model3b, VAS_Model4a)

Analysis of Variance Table

Model 1: VAS ~ Hull_bVSD_75 + Hull_b
Model 2: VAS ~ Hull_bVSD_75 + autoVSA + Hull_b
  Res.Df   RSS Df Sum of Sq      F  Pr(>F)  
1     37 25145                              
2     36 21202  1    3942.4 6.6939 0.01386 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Model 4b

## Specifying Model 4

VAS_Model4b <- lm(VAS ~ Hull_bVSD_75 + autoVSA, data = AcousticData)

## Model 4 Assumption Check

performance::check_model(VAS_Model4b)

performance::check_collinearity(VAS_Model4b)

# Check for Multicollinearity

Low Correlation

         Term  VIF Increased SE Tolerance
 Hull_bVSD_75 1.11         1.05      0.90
      autoVSA 1.11         1.05      0.90

## Model 4 Summary

summary(VAS_Model4b)

Call:
lm(formula = VAS ~ Hull_bVSD_75 + autoVSA, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-47.697 -17.980   9.493  20.999  36.709 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)   48.3844    14.1490   3.420  0.00154 **
Hull_bVSD_75   1.0096     2.0065   0.503  0.61783   
autoVSA        0.1236     1.3646   0.091  0.92834   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 26.85 on 37 degrees of freedom
Multiple R-squared:  0.00862,   Adjusted R-squared:  -0.04497 
F-statistic: 0.1608 on 2 and 37 DF,  p-value: 0.852

Model 5

## Specifying Model 5

VAS_Model5 <- lm(VAS ~ Hull_bVSD_75 + autoVSA + VSA_b, data = AcousticData)

## Model 5 Assumption Check

performance::check_model(VAS_Model5)

## Model 5 Summary

summary(VAS_Model5)

Call:
lm(formula = VAS ~ Hull_bVSD_75 + autoVSA + VSA_b, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-44.376 -19.256   5.331  18.101  42.048 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)   33.1721    13.2237   2.509  0.01677 * 
Hull_bVSD_75  -0.2154     1.8022  -0.120  0.90554   
autoVSA       -1.0281     1.2475  -0.824  0.41529   
VSA_b          6.7794     1.9786   3.426  0.00154 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 23.64 on 36 degrees of freedom
Multiple R-squared:  0.2524,    Adjusted R-squared:  0.1901 
F-statistic: 4.052 on 3 and 36 DF,  p-value: 0.01401

## Model 4 and Model 5 Comparison

anova(VAS_Model4b, VAS_Model5)

Analysis of Variance Table

Model 1: VAS ~ Hull_bVSD_75 + autoVSA
Model 2: VAS ~ Hull_bVSD_75 + autoVSA + VSA_b
  Res.Df   RSS Df Sum of Sq     F   Pr(>F)   
1     37 26675                               
2     36 20115  1    6559.9 11.74 0.001545 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Model 6

## Specifying Model 6

VAS_Model6 <- lm(VAS ~ Hull_bVSD_75 + autoVSA + VSA_b + vowel_ED_b, data = AcousticData)

## Model 6 Assumption Check

performance::check_model(VAS_Model6)

## Model 6 Summary

summary(VAS_Model6)

Call:
lm(formula = VAS ~ Hull_bVSD_75 + autoVSA + VSA_b + vowel_ED_b, 
    data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-42.388 -15.904   5.524  18.013  40.040 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept)  19.74662   22.78983   0.866    0.392  
Hull_bVSD_75  0.04266    1.84869   0.023    0.982  
autoVSA      -1.31095    1.31483  -0.997    0.326  
VSA_b         5.36014    2.79118   1.920    0.063 .
vowel_ED_b   10.79666   14.87603   0.726    0.473  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 23.8 on 35 degrees of freedom
Multiple R-squared:  0.2635,    Adjusted R-squared:  0.1793 
F-statistic: 3.131 on 4 and 35 DF,  p-value: 0.02658

## Model 5 and Model 6 Comparison

anova(VAS_Model5, VAS_Model6)

Analysis of Variance Table

Model 1: VAS ~ Hull_bVSD_75 + autoVSA + VSA_b
Model 2: VAS ~ Hull_bVSD_75 + autoVSA + VSA_b + vowel_ED_b
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1     36 20115                           
2     35 19817  1    298.25 0.5268 0.4728

Final Model

## Specifying Final Model

VAS_Model_final <- lm(VAS ~ VSA_b, data = AcousticData)

## Final Model Assumption Check

performance::check_model(VAS_Model_final)

## Final Model Summary

summary(VAS_Model_final)

Call:
lm(formula = VAS ~ VSA_b, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-44.956 -15.943   6.754  17.153  43.062 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   24.703      8.761   2.820  0.00760 **
VSA_b          6.163      1.798   3.427  0.00148 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 23.26 on 38 degrees of freedom
Multiple R-squared:  0.2361,    Adjusted R-squared:  0.216 
F-statistic: 11.74 on 1 and 38 DF,  p-value: 0.001482

confint(VAS_Model_final)

               2.5 %    97.5 %
(Intercept) 6.966936 42.438084
VSA_b       2.521887  9.803467

Research Q2: Relationship between OT and VAS

Model 1

# Specify Model

OT_VAS_model <- lm(transAcc ~ VAS*Etiology + VAS*Sex, data = AcousticData)

# Assumption Check

performance::check_model(OT_VAS_model)

# Model Results

summary(OT_VAS_model)

Final Linear Model

# Specify Final Model

OT_VAS_final <- lm(transAcc ~ VAS, data = AcousticData)

confint(OT_VAS_final)

# Model Results

summary(OT_VAS_final)

Corner Dispersion

Looking at corner dispersion as the sole predictor.

# Specify Final Model

OT_cornDisp <- lm(transAcc ~ vowel_ED_b, data = AcousticData)
summary(OT_cornDisp)

Call:
lm(formula = transAcc ~ vowel_ED_b, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-53.949 -10.348   4.268  14.982  25.903 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)    6.227     18.611   0.335  0.73977   
vowel_ED_b    25.564      8.946   2.857  0.00689 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21.98 on 38 degrees of freedom
Multiple R-squared:  0.1769,    Adjusted R-squared:  0.1552 
F-statistic: 8.165 on 1 and 38 DF,  p-value: 0.006891

VAS_cornDisp <- lm(VAS ~ vowel_ED_b, data = AcousticData)
summary(VAS_cornDisp)

Call:
lm(formula = VAS ~ vowel_ED_b, data = AcousticData)

Residuals:
    Min      1Q  Median      3Q     Max 
-52.146 -13.413   7.142  18.719  32.964 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   -2.859     20.636  -0.139   0.8905  
vowel_ED_b    26.819      9.920   2.704   0.0102 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 24.37 on 38 degrees of freedom
Multiple R-squared:  0.1613,    Adjusted R-squared:  0.1393 
F-statistic:  7.31 on 1 and 38 DF,  p-value: 0.0102

Manuscript Tables

##Descriptives Table

gtData <- AcousticData %>%
  rbind(.,AcousticData %>%
          dplyr::mutate(Etiology = "All Etiologies")) %>%
  rbind(.,AcousticData %>%
          rbind(.,AcousticData %>%
          dplyr::mutate(Etiology = "All Etiologies")) %>%
          dplyr::mutate(Sex = "All")) %>%
  dplyr::mutate(Sex = as.factor(Sex),
                Etiology = as.factor(Etiology)) %>%
  dplyr::group_by(Sex, Etiology) %>%
  dplyr::summarize(VSA_mean = mean(VSA_b, na.rm =T), VSA_sd = sd(VSA_b, na.rm = T),
                   Disp_mean = mean(vowel_ED_b, na.rm =T), Disp_sd = sd(vowel_ED_b, na.rm =T),
                   autoVSA_mean = mean(autoVSA, na.rm =T), autoVSA_sd = sd(autoVSA, na.rm = T),
                   Hull_mean = mean(Hull_b, na.rm =T), Hull_sd = sd(Hull_b, na.rm =T),
                   VSD25_mean = mean(Hull_bVSD_25, na.rm =T), VSD25_sd = sd(Hull_bVSD_25, na.rm =T),
                   VSD50_mean = mean(Hull_bVSD_50, na.rm =T), VSD50_sd = sd(Hull_bVSD_50, na.rm =T),
                   VSD75_mean = mean(Hull_bVSD_75, na.rm =T), VSD75_sd = sd(Hull_bVSD_75, na.rm =T),
                   VAS_mean = mean(VAS, na.rm =T), VAS_sd = sd(VAS, na.rm =T),
                   OT_mean = mean(transAcc, na.rm =T), OT_sd = sd(transAcc, na.rm =T)) %>%
  pivot_longer(cols = VSA_mean:OT_sd, names_to = "Measure",
               values_to = "Value") %>%
  dplyr::mutate(Value = round(Value, digits = 2),
                meanSD = ifelse(grepl("_mean",Measure),"M","sd"),
                Measure = gsub("_mean","",Measure),
                Measure = gsub("_sd","",Measure),
                Etiology = paste(Etiology,meanSD, sep = "_"),
                Sex = case_when(
                  Sex == "All" ~ "All Speakers",
                  Sex == "M" ~ "Male Speakers",
                  Sex == "F" ~ "Female Speakers"
                )) %>%
  dplyr::select(!meanSD) %>%
  pivot_wider(names_from = Etiology, values_from = "Value") %>%
  dplyr::filter(Measure != "VSD50")

gtData %>%
  gt::gt(
    rowname_col = "Measure",
    groupname_col = "Sex",
  ) %>%
  fmt_number(
    columns = 'All Etiologies_M':PD_sd,
    decimals = 2
  ) %>%
  tab_spanner(
    label = "All Etiologies",
    columns = c('All Etiologies_M', 'All Etiologies_sd')
  ) %>%
    tab_spanner(
    label = "ALS",
    columns = c(ALS_M, ALS_sd)
  ) %>%
  tab_spanner(
    label = "PD",
    columns = c(PD_M, PD_sd)
  ) %>%
  tab_spanner(
    label = "HD",
    columns = c(HD_M, HD_sd)
  ) %>%
  tab_spanner(
    label = "Ataxic",
    columns = c(Ataxic_M, Ataxic_sd)
  ) %>%
  gt::cols_move_to_start(
    columns = c('All Etiologies_M','All Etiologies_sd')
  ) %>%
  row_group_order(
    groups = c("All Speakers", "Female Speakers", "Male Speakers")
    ) %>%
  cols_label(
     'All Etiologies_M' = "M",
     'All Etiologies_sd' = "SD",
     ALS_M = "M",
     ALS_sd = "SD",
     PD_M = "M",
     PD_sd = "SD",
     HD_M = "M",
     HD_sd = "SD",
     Ataxic_M = "M",
     Ataxic_sd = "SD"
  ) %>%
  gtsave("DescriptivesTable.html", path = "Tables")

OT Model

sjPlot::tab_model(OT_Model1,OT_Model2,OT_Model3,
                  OT_Model4, OT_Model5,OT_Model6, OT_Model_final,
                  show.ci = F,
                  p.style = "stars",
                  file = "Tables/OT Models.html")

VAS Model

sjPlot::tab_model(VAS_Model1,
                  VAS_Model2,
                  VAS_Model3,
                  VAS_Model4,
                  VAS_Model5,
                  VAS_Model6,
                  VAS_Model_final,
                  show.ci = F,
                  p.style = "stars",
                  file = "Tables/VAS Models.html")

OT vs. VAS

sjPlot::tab_model(OT_VAS_model,OT_VAS_final,
                  show.ci = F,
                  show.reflvl = TRUE,
                  p.style = "stars",
                  file = "Tables/OT and VAS Comparison.html")

Manuscript Figures

Example Measures

formantColor <- "grey"
formantAlpha <- .95
lineColor <- "white"
lineAlpha <- .8

vowelData <- rio::import("Prepped Data/Vowel Data.csv") %>%
  dplyr::filter(Speaker == "AF8")

  Pitch_PRAAT <- list.files(path = paste("Prepped Data/Example Data/", sep = ""), 
                              pattern = ".Pitch", ignore.case = T) %>%
    paste("Prepped Data/Example Data/",., sep = "") %>%
    read.delim(., header = F) %>%
    dplyr::rename(Pitch = V1) %>%
    dplyr::mutate(Pitch = gsub("--undefined--",NA,Pitch),
                  Pitch = as.numeric(Pitch))
  
  Formants_PRAAT <- list.files(path = paste("Prepped Data/Example Data/", sep = ""), 
                              pattern = "_Formant", ignore.case = T) %>%
    paste("Prepped Data/Example Data/",., sep = "") %>%
    read.delim(., header = T) %>%
    dplyr::select(!c(nformants, B1.Hz., B2.Hz., B3.Hz., F4.Hz., B4.Hz., F5.Hz., B5.Hz.)) %>%
    dplyr::rename(Time_s = time.s.,
                  F1_Hz = F1.Hz.,
                  F2_Hz = F2.Hz.,
                  F3_Hz = F3.Hz.) %>%
    dplyr::mutate(F1_Hz = ifelse(F1_Hz == 0, NA, F1_Hz),
                  F2_Hz = ifelse(F2_Hz == 0, NA, F2_Hz),
                  F3_Hz = ifelse(F3_Hz == 0, NA, F3_Hz)) %>%
    dplyr::mutate(F1_Hz = as.numeric(F1_Hz),
                  F2_Hz = as.numeric(F2_Hz),
                  F3_Hz = suppressWarnings(as.numeric(F3_Hz)),
                  Time_ms = Time_s / 1000,
                  F1_kHz = F1_Hz / 1000,
                  F2_kHz = F2_Hz / 1000,
                  F3_kHz = F3_Hz / 1000) %>%
    dplyr::select(!Time_s) %>%
    dplyr::relocate(Time_ms, .before = F1_Hz) %>%
    cbind(.,Pitch_PRAAT)
  
  c <- 2
  while(c < NROW(Formants_PRAAT)){
    Formants_PRAAT$F1_Hz[c] <- ifelse(is.na(Formants_PRAAT$F1_Hz[c-1]) &&
                                        is.na(Formants_PRAAT$F1_Hz[c+1]),
                                      NA,
                                      Formants_PRAAT$F1_Hz[c])
    Formants_PRAAT$F2_Hz[c] <- ifelse(is.na(Formants_PRAAT$F2_Hz[c-1]) &&
                                        is.na(Formants_PRAAT$F2_Hz[c+1]),
                                      NA,
                                      Formants_PRAAT$F2_Hz[c])
    c <- c + 1
  }
  rm(c)
  
  Formants_PRAAT <- Formants_PRAAT %>%
    dplyr::filter(!is.na(Pitch)) %>%
    dplyr::mutate(F1_mad = (abs(F1_Hz - median(F1_Hz))/ mad(F1_Hz, constant = 1.4826)) > 2.5,
                  F2_mad = (abs(F2_Hz - median(F2_Hz))/ mad(F2_Hz, constant = 1.4826)) > 2.5) %>%
    dplyr::filter(F1_mad == FALSE & F2_mad == FALSE) %>%
    dplyr::mutate(mDist = mahalanobis(cbind(.$F1_Hz, .$F2_Hz),
                                      colMeans(cbind(.$F1_Hz, .$F2_Hz)),
                                      cov = cov(cbind(.$F1_Hz, .$F2_Hz))),
                  mDist_sd = abs(scale(mDist,center = T))) %>%
    dplyr::filter(mDist_sd < 2) %>%
    dplyr::select(!c(F1_mad, F2_mad, mDist, mDist_sd)) %>%
    dplyr::mutate(F1_z = scale(F1_Hz, center = TRUE),
                  F2_z = scale(F2_Hz, center = TRUE),
                  F3_z = scale(F3_Hz, center = TRUE),
                  F1_b = emuR::bark(F1_Hz),
                  F2_b = emuR::bark(F2_Hz),
                  F3_b = emuR::bark(F3_Hz))
  
  rm(Pitch_PRAAT)
  
  
## Corner Dispersion ----
  wedge <- vowelData %>%
    dplyr::group_by(Vowel) %>%
    dplyr::summarize(mean_F1 = mean(F1_tempMid),
              mean_F2 = mean(F2_tempMid),
              mean_F1_z = mean(F1_z_tempMid),
              mean_F2_z = mean(F2_z_tempMid),
              mean_F1_b = mean(F1_b_tempMid),
              mean_F2_b = mean(F2_b_tempMid)) %>%
    dplyr::filter(Vowel == "v")
    
  corner_dis <- vowelData %>%
    dplyr::filter(Vowel != "v") %>%
    dplyr::group_by(Vowel) %>%
    dplyr::summarize(mean_F1 = mean(F1_tempMid),
              mean_F2 = mean(F2_tempMid),
              mean_F1_z = mean(F1_z_tempMid),
              mean_F2_z = mean(F2_z_tempMid),
              mean_F1_b = mean(F1_b_tempMid),
              mean_F2_b = mean(F2_b_tempMid)) %>%
    dplyr::mutate(Vowel_ED = sqrt((mean_F1-wedge$mean_F1)^2 + (mean_F2-wedge$mean_F2)^2),
                  Vowel_ED_z = sqrt((mean_F1_z-wedge$mean_F1_z)^2 + (mean_F2_z-wedge$mean_F2_z)^2),
                  Vowel_ED_b = sqrt((mean_F1_b-wedge$mean_F1_b)^2 + (mean_F2_b-wedge$mean_F2_b)^2))

    
# Plot Corner Dispersion
      # Changing to IPA symbols
      corner_dis <- corner_dis %>%
        dplyr::mutate(Vowel = dplyr::case_when(
          Vowel == "ae" ~ "æ",
          TRUE ~ Vowel
        ))
      
      wedge <- wedge %>%
        dplyr::mutate(Vowel = case_when(
          Vowel == "v" ~ "ʌ",
          TRUE ~ Vowel
        ))
      
      CDplot <- ggplot(aes(x=F2_b,
                           y=F1_b),
                       data = Formants_PRAAT,
                       inherit.aes = FALSE) + 
      geom_point(shape = 21,
                 alpha = formantAlpha,
                 color = formantColor) + 
      geom_line(aes(x = mean_F2_b,
                    y = mean_F1_b),
                data = corner_dis %>%
                  dplyr::select(Vowel:mean_F2_b) %>%
                  dplyr::filter(Vowel == "i") %>%
                  rbind(.,wedge),
                color = lineColor,
                size = 1.5,
                alpha = lineAlpha) +
      geom_line(aes(x = mean_F2_b,
                    y = mean_F1_b),
                data = corner_dis %>%
                  dplyr::select(Vowel:mean_F2_b) %>%
                  dplyr::filter(Vowel == "a") %>%
                  rbind(.,wedge),
                color = lineColor,
                size = 1.5,
                alpha = lineAlpha) +
      geom_line(aes(x = mean_F2_b,
                    y = mean_F1_b),
                data = corner_dis %>%
                  dplyr::select(Vowel:mean_F2_b) %>%
                  dplyr::filter(Vowel == "æ") %>%
                  rbind(.,wedge),
                color = lineColor,
                size = 1.5,
                alpha = lineAlpha) +
      geom_line(aes(x = mean_F2_b,
                    y = mean_F1_b),
                data = corner_dis %>%
                  dplyr::select(Vowel:mean_F2_b) %>%
                  dplyr::filter(Vowel == "u") %>%
                  rbind(.,wedge),
                color = lineColor,
                size = 1.5,
                alpha = lineAlpha) +
      geom_point(aes(x = mean_F2_b,
                     y = mean_F1_b,
                     color = Vowel),
                 data = corner_dis %>%
                  dplyr::select(Vowel:mean_F2_b) %>%
                  rbind(.,wedge),
                 inherit.aes = FALSE,
                 size = 5) +
      scale_y_reverse() +
      scale_x_reverse() +
      theme_classic() + labs(title = paste("Corner Dispersion")) + xlab("F2 (bark)") + ylab("F1 (bark)") +
      theme(plot.title = element_text(hjust = 0.5),
            aspect.ratio = 1) +
          scale_color_manual(values = c("a" = "#1AAD77",
                                        "æ" = "#1279B5",
                                        "i" = "#FFBF00",
                                        "u" = "#FD7853",
                                        "ʌ" = "#BF3178"))
    CDplot
    
      rm(corner_dis, wedge)
      
## Vowel Space Area ----
  VSA_coords <- vowelData %>%
    dplyr::filter(Vowel != "v") %>%
    dplyr::group_by(Vowel) %>%
    dplyr::summarize(mean_F1 = mean(F1_tempMid),
              mean_F2 = mean(F2_tempMid),
              mean_F1_z = mean(F1_z_tempMid),
              mean_F2_z = mean(F2_z_tempMid),
              mean_F1_b = mean(F1_b_tempMid),
              mean_F2_b = mean(F2_b_tempMid)) 
  
### Plotting VSA
    VSA_coords <- VSA_coords %>%
        dplyr::mutate(Vowel = case_when(
          Vowel == "ae" ~ "æ",
          TRUE ~ Vowel
        ))
    
    VSAplot <- ggplot(aes(x = F2_b,
                          y = F1_b),
                      data = Formants_PRAAT,
                      inherit.aes = FALSE) + 
      geom_point(shape = 21,
                 alpha = formantAlpha,
                 color = formantColor) + 
      geom_polygon(aes(x = mean_F2_b,
                       y = mean_F1_b),
                   data = VSA_coords,
                   alpha = lineAlpha,
                   color = lineColor,
                   fill=NA,
                   size = 1.5) +
      geom_point(aes(x = mean_F2_b,
                     y = mean_F1_b,
                     color = Vowel),
                 data = VSA_coords,
                 inherit.aes = FALSE,
                 size = 5) +
      scale_y_reverse() +
      scale_x_reverse() +
      guides(color = FALSE) +
      theme_classic() + labs(title = "VSA") + xlab("F2 (bark)") + ylab("F1 (bark)") +
      theme(plot.title = element_text(hjust = 0.5),
            aspect.ratio = 1) +
                scale_color_manual(values = c("a" = "#1AAD77",
                                        "æ" = "#1279B5",
                                        "i" = "#FFBF00",
                                        "u" = "#FD7853"))
    VSAplot
  
  rm(VSA_coords)
  
## Automatic VSA ----
  RefData <- openxlsx::read.xlsx("Prepped Data/Example Data/Hillenbrand Vowel Data.xlsx",
                     sheet = "Hillenbrand Vowel Data") %>%
  dplyr::rename(Speaker = 1,
                F1_Hz = 4,
                F2_Hz = 5) %>%
  dplyr::select(Speaker, F1_Hz, F2_Hz) %>%
  dplyr::filter(!grepl("b|g",Speaker)) %>%
  dplyr::mutate(F1 = emuR::bark(F1_Hz),
                F2 = emuR::bark(F2_Hz),
                Sex = ifelse(grepl("m",Speaker),"M","F"),
                Vowel = case_when(
                  grepl(pattern = "uw", Speaker) ~ "u",
                  grepl(pattern = "ah", Speaker) ~ "a",
                  grepl(pattern = "iy", Speaker) ~ "i",
                  grepl(pattern = "uh", Speaker) ~ "v",
                  grepl(pattern = "ae", Speaker) ~ "ae",
                  grepl(pattern = "aw", Speaker) ~ "ɔ",
                  grepl(pattern = "eh", Speaker) ~ "ɛ",
                  grepl(pattern = "er", Speaker) ~ "ɝ",
                  grepl(pattern = "ei", Speaker) ~ "eɪ",
                  grepl(pattern = "ih", Speaker) ~ "ɪ",
                  grepl(pattern = "oa", Speaker) ~ "o",
                  grepl(pattern = "oo", Speaker) ~ "ʊ"
                )) %>%
  dplyr::filter(!is.na(Vowel)) %>%
  dplyr::select(!Speaker) %>%
  dplyr::group_by(Sex, Vowel) %>%
  dplyr::summarise(F1_Ref = mean(F1),
                   F2_Ref = mean(F2)) %>%
    dplyr::ungroup() %>%
    dplyr::filter(Sex == "F") %>%
    dplyr::select(!c(Sex, Vowel))
  
  autoVSA <- Formants_PRAAT %>%
    dplyr::select(F1_b,F2_b) %>%
    stats::kmeans(., centers = RefData, nstart = 100, iter.max = 1000) %>%
    .$centers %>%
    as.data.frame()
  
  convexCoords <- autoVSA %>%
    as.matrix() %>%
    grDevices::chull()
  
  autoConvex <- autoVSA %>%
        slice(convexCoords)
  
### Plotting Automatic VSA
    autoVSAplot <- ggplot(aes(x = F2_b,
                          y = F1_b),
                      data = Formants_PRAAT,
                      inherit.aes = FALSE) + 
      geom_point(shape = 21,
                 alpha = formantAlpha,
                 color = formantColor) + 
      geom_polygon(aes(x = F2_b,
                       y = F1_b),
                   data = autoConvex,
                   alpha = lineAlpha,
                   color = lineColor,
                   fill=NA,
                   size = 1.5) +
      geom_point(aes(x = F2_b,
                     y = F1_b),
                 data = autoVSA,
                 inherit.aes = FALSE,
                 size = 2) +
      scale_y_reverse() +
      scale_x_reverse() +
      #guides(color = FALSE) +
      theme_classic() +
      labs(title = "Automatic VSA") +
      xlab("F2 (bark)") +
      ylab("F1 (bark)") +
      theme(plot.title = element_text(hjust = 0.5),
            aspect.ratio = 1)
    autoVSAplot
  
  rm(VSA_coords)
  
## Hull ----
### Plotting Hull
      convexCoords <- Formants_PRAAT %>%
        dplyr::select(F1_b, F2_b) %>%
        as.matrix() %>%
        grDevices::chull()
      convex <- Formants_PRAAT %>%
        slice(convexCoords)

      hullPlot <- ggplot(aes(F2_b, F1_b),
                         data = Formants_PRAAT) +
        geom_point(shape = 21,
                 alpha = formantAlpha,
                 color = formantColor) +
        geom_polygon(data = convex,
                     alpha = .5,
                     color = "#1279B5",
                     fill = NA,
                     size = 1.5) +
        scale_y_reverse() +
        scale_x_reverse() +
        theme_classic() + labs(title = "VSA Hull") + xlab("F2 (bark)") + ylab("F1 (bark)") +
      theme(plot.title = element_text(hjust = 0.5),
            aspect.ratio = 1)
      hullPlot
    
  
## Vowel Space Density ----

# Bark Normalized Density ----
# selecting the bandwidth
H_hpi <- ks::Hpi(x = Formants_PRAAT[,c("F2_b","F1_b")], pilot = "samse", pre = "scale", binned = T)

# compute 2d kde
k <- kde(x = Formants_PRAAT[,c("F2_b","F1_b")],
         H = H_hpi,
         binned = T,
         gridsize = 250)

#density <- k[["estimate"]]

# Before we can plot the density estimate we need to melt it into long format
mat.melted <- data.table::melt(k$estimate)
names(mat.melted) <- c("x", "y", "density")

# We need to add two more colums to preserve the axes units
mat.melted$F2.b <- rep(k$eval.points[[1]], times = nrow(k$estimate))
mat.melted$F1.b <- rep(k$eval.points[[2]], each = nrow(k$estimate))
mat.melted$density <- scales::rescale(mat.melted$density, to = c(0, 1))

# VSD - 25
nVSD_25 <- mat.melted %>%
  dplyr::filter(density > .25) %>%
  dplyr::select(F2.b,F1.b, density) %>%
  dplyr::rename(Density = density)

convexCoords <- nVSD_25 %>%
  dplyr::select(F2.b, F1.b) %>%
  as.matrix() %>%
  #grDevices::xy.coords() %>%
  grDevices::chull()
nconvex_25 <- nVSD_25 %>%
  slice(convexCoords)

# VSD - 75
nVSD_75 <- mat.melted %>%
  dplyr::filter(density > .75) %>%
  dplyr::select(F2.b,F1.b, density) %>%
  dplyr::rename(Density = density)

convexCoords <- nVSD_75 %>%
  dplyr::select(F2.b, F1.b) %>%
  as.matrix() %>%
  grDevices::chull()
nconvex_75 <- nVSD_75 %>%
  slice(convexCoords)

# Plotting Z Normalized VSD 
    rf <- colorRampPalette(rev(RColorBrewer::brewer.pal(11, "Spectral")))
    r <- rf(32)
    
    plotData <- mat.melted %>%
                        dplyr::rename(Density = density) %>%
                        dplyr::mutate(VSDlabel = dplyr::case_when(
                          Density < .25 ~ "none",
                          Density > .25 && Density < .75 ~ "VSD25",
                          TRUE ~ "VSD75"
                        ))
geom.text.size <- 2
    VSDplot <- ggplot(data = plotData,
                      aes(x = F2.b,
                          y = F1.b,
                          fill = Density)) + 
      geom_tile() + 
      scale_fill_viridis_c() +
      scale_x_reverse(expand = c(0, 0), 
                      breaks = round(seq(min(mat.melted$F2.b), 
                                         max(mat.melted$F2.b), by = 2))) +
      scale_y_reverse(expand = c(0, 0),
                      breaks = round(seq(min(mat.melted$F1.b),
                                         max(mat.melted$F1.b), by = 2))) + 
      ylab("F1 (bark)") + xlab("F2 (bark)") +
      labs(title = "VSD Hull") +
      theme(plot.title = element_text(hjust = 0.5),
            aspect.ratio = 1) +
      geom_polygon(data = nconvex_25, alpha = lineAlpha, color = lineColor, size = 1.5, fill = NA, linetype = 2) +
      geom_polygon(data = nconvex_75, alpha = lineAlpha, color = lineColor, size = 1.5, fill = NA, linetype = 1) +
    # VSD 25 Label
      annotate(geom = "curve",
               x = 6.9, y = 1.7+.5,
               xend = 8.5, yend = 3.5,
               curvature = -.3,
               arrow = arrow(length = unit(2, "mm")),
               color = "white") +
      annotate(geom = "text",
               x = 7.5, y = 1.7,
               label = deparse(bquote(VSD[25])),
               hjust = "center",
               color = "white",
               parse=TRUE) +
    # VSD 75 Label
      annotate(geom = "curve",
               x = 7.5, y = 7.5-.5,
               xend = 11.35, yend = 5.5,
               curvature = .3,
               arrow = arrow(length = unit(2, "mm")),
               color = "white") +
      annotate(geom = "text",
               x = 7, y = 7.5,
               label = deparse(bquote(VSD[75])),
               hjust = "center",
               color = "white",
               parse = TRUE)
     VSDplot

# Combined Plot
     
     row1 <- VSAplot + CDplot +
        patchwork::plot_layout(guides = 'collect',
                         ncol = 2) & theme(legend.position = 'right')
     row2 <- autoVSAplot + hullPlot + VSDplot +
        patchwork::plot_layout(guides = 'collect',
                         ncol = 3) & theme(legend.position = 'right')
     
     measuresPlot <- row1 / row2 + patchwork::plot_layout(heights = c(1/2, 1/2), byrow = FALSE)
     measuresPlot
     
     rm(row1, row2)

ggsave(filename = "Plots/Measures.png",
       plot = measuresPlot,
       height = 5.5,
       width = 8,
       scale = .8)

Filtering Process

formantAlpha <- .20
myPal <- c("#1279B5","#2D2D37")

Pitch_PRAAT <- list.files(path = paste("Prepped Data/Example Data/", sep = ""), 
                              pattern = ".Pitch", ignore.case = T) %>%
    paste("Prepped Data/Example Data/",., sep = "") %>%
    read.delim(., header = F) %>%
    dplyr::rename(Pitch = V1) %>%
    dplyr::mutate(Pitch = gsub("--undefined--",NA,Pitch),
                  Pitch = as.numeric(Pitch))

Formants_PRAAT <- list.files(path = paste("Prepped Data/Example Data/", sep = ""), 
                              pattern = "_Formant", ignore.case = T) %>%
    paste("Prepped Data/Example Data/",., sep = "") %>%
    read.delim(., header = T) %>%
    dplyr::select(!c(nformants, B1.Hz., B2.Hz., B3.Hz., F4.Hz., B4.Hz., F5.Hz., B5.Hz.)) %>%
    dplyr::rename(Time_s = time.s.,
                  F1_Hz = F1.Hz.,
                  F2_Hz = F2.Hz.,
                  F3_Hz = F3.Hz.) %>%
    dplyr::mutate(F1_Hz = ifelse(F1_Hz == 0, NA, F1_Hz),
                  F2_Hz = ifelse(F2_Hz == 0, NA, F2_Hz),
                  F3_Hz = ifelse(F3_Hz == 0, NA, F3_Hz)) %>%
    dplyr::mutate(F1_Hz = as.numeric(F1_Hz),
                  F2_Hz = as.numeric(F2_Hz),
                  F3_Hz = suppressWarnings(as.numeric(F3_Hz)),
                  Time_ms = Time_s / 1000,
                  F1_kHz = F1_Hz / 1000,
                  F2_kHz = F2_Hz / 1000,
                  F3_kHz = F3_Hz / 1000,
                  F1_b = emuR::bark(F1_Hz),
                  F2_b = emuR::bark(F2_Hz),
                  F3_b = emuR::bark(F3_Hz)) %>%
    dplyr::select(!Time_s) %>%
    dplyr::relocate(Time_ms, .before = F1_Hz) %>%
    cbind(.,Pitch_PRAAT)
  
  c <- 2
  while(c < NROW(Formants_PRAAT)){
    Formants_PRAAT$F1_Hz[c] <- ifelse(is.na(Formants_PRAAT$F1_Hz[c-1]) &&
                                        is.na(Formants_PRAAT$F1_Hz[c+1]),
                                      NA,
                                      Formants_PRAAT$F1_Hz[c])
    Formants_PRAAT$F2_Hz[c] <- ifelse(is.na(Formants_PRAAT$F2_Hz[c-1]) &&
                                        is.na(Formants_PRAAT$F2_Hz[c+1]),
                                      NA,
                                      Formants_PRAAT$F2_Hz[c])
    c <- c + 1
  }
  rm(c)
  
  # Raw Formants ----
  f1 <- ggplot(aes(x=F2_b,
                   y=F1_b),
               data = Formants_PRAAT) + 
      geom_point(shape = 21, color = myPal[2]) +
      scale_y_reverse() +
      scale_x_reverse() +
      scale_color_manual(values = myPal) +
      theme_classic() + labs(title = paste("Raw Formant Values")) + xlab("F2 (bark)") + ylab("F1 (bark)") +
      theme(plot.title = element_text(hjust = 0.5),
            aspect.ratio = 1,
            legend.title = element_blank(),
            legend.text = element_text(size=12))
    
# Step #1: Voiced Segments ----
    plotData <- Formants_PRAAT %>%
                   dplyr::mutate(isOutlier = case_when(
                     is.na(Pitch) ~ "Removed",
                     TRUE ~ "Retained"
                   ))
    f2 <- ggplot(data = plotData,
                 aes(x = F2_b,
                     y = F1_b,
                     color = isOutlier)) + 
      geom_point(shape = 21, data = plotData %>%
                   dplyr::filter(isOutlier == "Removed")) +
      geom_point(shape = 21, data = plotData %>%
             dplyr::filter(isOutlier == "Retained")) +
      scale_y_reverse() +
      scale_x_reverse() +
      scale_color_manual(values = myPal) +
      theme_classic() + labs(title = paste("Step #1:\nVoiced Segments")) +
      xlab("F2 (bark)") +
      ylab("F1 (bark)") +
      theme(plot.title = element_text(hjust = 0.5),
            aspect.ratio = 1,
            legend.title = element_blank(),
            legend.text = element_text(size=12))
    
# Step 2: MAD ----
    plotData <- Formants_PRAAT %>%
      dplyr::filter(!is.na(Pitch)) %>%
      dplyr::mutate(F1_mad = (abs(F1_Hz - median(F1_Hz))/ mad(F1_Hz, constant = 1.4826)) > 2.5,
                    F2_mad = (abs(F2_Hz - median(F2_Hz))/ mad(F2_Hz, constant = 1.4826)) > 2.5,
                    isOutlier = case_when(
                      F1_mad == TRUE | F2_mad == TRUE ~ "Removed",
                      TRUE ~ "Retained"
               ))
    
    f3 <- ggplot(data = plotData,
                 aes(x = F2_b,
                     y = F1_b,
                     color = isOutlier)) + 
      geom_point(shape = 21, data = plotData %>%
                   dplyr::filter(isOutlier == "Removed")) +
      geom_point(shape = 21, data = plotData %>%
             dplyr::filter(isOutlier == "Retained")) +
      scale_y_reverse() +
      scale_x_reverse() +
      scale_color_manual(values = myPal) +
      theme_classic() +
      labs(title = paste("Step #2:\nMedian Absolute Deviation")) +
      xlab("F2 (bark)") +
      ylab("F1 (bark)") +
      theme(plot.title = element_text(hjust = 0.5),
            aspect.ratio = 1,
            legend.title = element_blank(),
            legend.text = element_text(size=12))
    
# Step 3: Mahalanhobis Distance ----
  plotData <- Formants_PRAAT %>%
      dplyr::filter(!is.na(Pitch)) %>%
      dplyr::mutate(F1_mad = (abs(F1_Hz - median(F1_Hz))/ mad(F1_Hz, constant = 1.4826)) > 2.5,
                    F2_mad = (abs(F2_Hz - median(F2_Hz))/ mad(F2_Hz, constant = 1.4826)) > 2.5) %>%
      dplyr::filter(F1_mad == FALSE & F2_mad == FALSE) %>%
      dplyr::mutate(mDist = mahalanobis(cbind(.$F1_Hz, .$F2_Hz),
                                      colMeans(cbind(.$F1_Hz, .$F2_Hz)),
                                      cov = cov(cbind(.$F1_Hz, .$F2_Hz))),
                  mDist_sd = abs(scale(mDist,center = T)),
                  isOutlier = case_when(
                    mDist_sd < 2 ~ "Retained",
                    TRUE ~ "Removed"
                  ))
    
    f4 <- ggplot(data = plotData,
                 aes(x = F2_b,
                     y = F1_b,
                     color = isOutlier)) + 
      geom_point(shape = 21, data = plotData %>%
                   dplyr::filter(isOutlier == "Removed")) +
      geom_point(shape = 21, data = plotData %>%
             dplyr::filter(isOutlier == "Retained")) +
      scale_y_reverse() +
      scale_x_reverse() +
      scale_color_manual(values = myPal) +
      theme_classic() + labs(title = paste("Step #3:\nMahalanobis Distance")) + xlab("F2 (bark)") + ylab("F1 (bark)") +
      theme(plot.title = element_text(hjust = 0.5),
            aspect.ratio = 1,
            legend.title = element_blank(),
            legend.text = element_text(size=12))
    
# Final Formants ----
    plotData <- Formants_PRAAT %>%
    dplyr::filter(!is.na(Pitch)) %>%
    dplyr::mutate(F1_mad = (abs(F1_Hz - median(F1_Hz))/ mad(F1_Hz, constant = 1.4826)) > 2.5,
                  F2_mad = (abs(F2_Hz - median(F2_Hz))/ mad(F2_Hz, constant = 1.4826)) > 2.5) %>%
    dplyr::filter(F1_mad == FALSE & F2_mad == FALSE) %>%
    dplyr::mutate(mDist = mahalanobis(cbind(.$F1_Hz, .$F2_Hz),
                                      colMeans(cbind(.$F1_Hz, .$F2_Hz)),
                                      cov = cov(cbind(.$F1_Hz, .$F2_Hz))),
                  mDist_sd = abs(scale(mDist,center = T))) %>%
    dplyr::filter(mDist_sd < 2)
    
    f5 <- ggplot(aes(x=F2_b,
                     y=F1_b),
                       data = plotData,
                       inherit.aes = FALSE) + 
      geom_point(shape = 21, color = myPal[2]) + 
      scale_y_reverse() +
      scale_x_reverse() +
      theme_classic() + labs(title = paste("Final Formant Values")) + xlab("F2 (bark)") + ylab("F1 (bark)") +
      theme(plot.title = element_text(hjust = 0.5),
            aspect.ratio = 1, legend.title = element_blank())
    
# Comibing plots
    filteredPlot <- f1 + f2 + f3 + f4 + f5 + patchwork::guide_area() +
      patchwork::plot_layout(guides = 'collect',
                         ncol = 3)
    filteredPlot
    
    ggsave(plot = filteredPlot, "Plots/Filtered Formants.png",
           height = 6,
           width = 8,
           units = "in",
           scale = .9)
  

OT vs. VAS

plotData_Int <- AcousticData %>%
  dplyr::filter(!grepl("_rel", Speaker)) %>%
  dplyr::group_by(Speaker) %>%
  dplyr::mutate(segMin = base::min(VAS, transAcc),
                segMax = base::max(VAS, transAcc),
                ratingAvg = mean(VAS, transAcc, na.rm = T),
                Speaker = as.factor(Speaker),
                Etiology = as.factor(Etiology)) %>%
  arrange(segMax)

my_pal <- c("#f26430", "#272D2D","#256eff")
# With a bit more style
plot_Int <- ggplot(plotData_Int) +
  geom_segment(aes(x = fct_inorder(Speaker),
                   xend = Speaker,
                   y = segMin,
                   yend = segMax,
                   color = Etiology)) +
  geom_point(aes(x = Speaker,
                 y = VAS,
                 color = Etiology),
             #color = my_pal[1],
             size = 3,
             shape = 19) +
  geom_point(aes(x = Speaker,
                 y = transAcc,
                 color = Etiology),
             #color = my_pal[2],
             size = 3,
             shape = 15) +
  coord_flip()+
  theme_classic() +
  theme(
    legend.position = "none",
    panel.border = element_blank(),
  ) +
  xlab("") +
  ylab("Speech Intelligibility") +
  ggtitle("Speech Intelligibility") +
  ylim(c(0,100))
plot_Int

myPal <- c("#1AAD77", "#1279B5", "#FFBF00", "#FD7853", "#BF3178")
myShapes <- c(16, 18, 17, 15)

scatter1 <- ggplot(plotData_Int,
                  aes(x = VAS,
                      y = transAcc,
                      color = Etiology,
                      shape = Etiology,
                      linetype = Etiology)) +
  geom_point() +
  #geom_smooth(method = "lm", se = F) +
  geom_abline(intercept = 0, slope = 1) +
  coord_cartesian(xlim = c(0,100), ylim = c(0,100)) +
  labs(x = "Intelligibility (VAS)", y = "Intelligibility (OT)") +
  scale_color_manual(values = myPal) +
  scale_shape_manual(values = myShapes) +
  theme_classic() +
  theme(aspect.ratio=1,
        legend.position="right")

scatter2 <- ggplot(plotData_Int,
                  aes(x = VAS,
                      y = transAcc,
                      color = Etiology,
                      shape = Etiology,
                      linetype = Etiology)) +
  #geom_point() +
  geom_smooth(method = "lm", se = F) +
  geom_abline(intercept = 0, slope = 1) +
  coord_cartesian(xlim = c(0,100), ylim = c(0,100)) +
  labs(x = "Intelligibility (VAS)", y = "Intelligibility (OT)") +
  scale_color_manual(values = myPal) +
  scale_shape_manual(values = myShapes) +
  theme_classic() +
  theme(aspect.ratio=1,
        legend.position="right")

combinedScatter <- ggarrange(scatter1,
                             scatter2,
                             common.legend = F,
                             ncol = 2,
                             nrow = 1)
combinedScatter

ggsave(filename = "Plots/OT and VAS Scatterplot.png",
       plot = combinedScatter,
       height = 2,
       width = 6,
       units = "in",
       scale = 1)

rm(scatter1, scatter2, combinedScatter)

Model Scatterplots

modelFigureData <- AcousticData %>%
  dplyr::filter(!grepl("_rel",Speaker)) %>%
  dplyr::select(Speaker, Etiology, Sex, autoVSA, VSA_b, vowel_ED_b, Hull_b, Hull_bVSD_25, Hull_bVSD_75, VAS, transAcc) %>%
  dplyr::mutate(Speaker = as.factor(Speaker),
                Etiology = as.factor(Etiology),
                Sex = as.factor(Sex)) %>%
  tidyr::pivot_longer(cols = VAS:transAcc, names_to = "IntType", values_to = "Int") %>%
  dplyr::mutate(IntType = case_when(
    IntType == "transAcc" ~ "OT",
    TRUE ~ "VAS"
  ),
                IntType = as.factor(IntType))

ylabel <- "Intelligibility"
myPal <- c("#2D2D37", "#1279B5")
myPalShape <- c(19, 1)

VSA <- modelFigureData %>%
  ggplot() +
  aes(x = VSA_b,
      y = Int,
      color = IntType,
      shape = IntType,
      linetype = IntType) +
  geom_point() +
  geom_smooth(method = "lm", se = T, fill = "light grey") +
  xlab(expression("VSA (Bark"^2*")")) +
  ylab(ylabel) +
  coord_cartesian(ylim = c(0,100)) +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        aspect.ratio=1) +
  scale_color_manual(values = myPal) +
  scale_shape_manual(values = myPalShape) +
  labs(color="Intelligibility Type",
       shape = "Intelligibility Type",
       linetype = "Intelligibility Type")

autoVSA <- modelFigureData %>%
  ggplot() +
  aes(x = autoVSA,
      y = Int,
      color = IntType,
      shape = IntType,
      linetype = IntType) +
  geom_point() +
  geom_smooth(method = "lm", se = T, fill = "light grey") +
  xlab(expression("VSA"[Automatic]*" (Bark"^2*")")) +
  ylab(ylabel) +
  coord_cartesian(ylim = c(0,100)) +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        aspect.ratio=1) +
  scale_color_manual(values = myPal) +
  scale_shape_manual(values = myPalShape) +
  labs(color="Intelligibility Type",
       shape = "Intelligibility Type",
       linetype = "Intelligibility Type")

disp <- modelFigureData %>%
  ggplot() +
  aes(x = vowel_ED_b,
      y = Int,
      color = IntType,
      shape = IntType,
      linetype = IntType) +
  geom_point() +
  geom_smooth(method = "lm", se = T, fill = "light grey") +
  xlab("Corner Dispersion (Bark)") +
  ylab(ylabel) +
  coord_cartesian(ylim = c(0,100)) +
  theme_classic() +
  theme(aspect.ratio=1) +
  scale_color_manual(values = myPal) +
  scale_shape_manual(values = myPalShape) +
  labs(color="Intelligibility Type",
       shape = "Intelligibility Type",
       linetype = "Intelligibility Type")

Hull <- modelFigureData %>%
  ggplot() +
  aes(x = Hull_b,
      y = Int,
      color = IntType,
      shape = IntType,
      linetype = IntType) +
  geom_point() +
  geom_smooth(method = "lm", se = T, fill = "light grey") +
  xlab(expression("VSA"[Hull]*" (Bark"^2*")")) +
  ylab(ylabel) +
  coord_cartesian(ylim = c(0,100)) +
  theme_classic() +
  theme(aspect.ratio=1) + theme(legend.position = "none") +
  scale_color_manual(values = myPal) +
  scale_shape_manual(values = myPalShape) +
  labs(color="Intelligibility Type",
       shape = "Intelligibility Type",
       linetype = "Intelligibility Type")

vsd25 <- modelFigureData %>%
  ggplot() +
  aes(x = Hull_bVSD_25,
      y = Int,
      color = IntType,
      shape = IntType,
      linetype = IntType) +
  geom_point() +
  geom_smooth(method = "lm", se = T, fill = "light grey") +
  xlab(expression("VSD"[25]*" (Bark"^2*")")) +
  ylab(ylabel) +
  coord_cartesian(ylim = c(0,100)) +
  theme_classic() +
  theme(aspect.ratio=1) + theme(legend.position = "none") +
  scale_color_manual(values = myPal) +
  scale_shape_manual(values = myPalShape) +
  labs(color="Intelligibility Type",
       shape = "Intelligibility Type",
       linetype = "Intelligibility Type")

vsd75 <- modelFigureData %>%
  ggplot() +
  aes(x = Hull_bVSD_75,
      y = Int,
      color = IntType,
      shape = IntType,
      linetype = IntType) +
  geom_point() +
  geom_smooth(method = "lm", se = T, fill = "light grey") +
  xlab(expression("VSD"[75]*" (Bark"^2*")")) +
  ylab(ylabel) +
  coord_cartesian(ylim = c(0,100)) +
  theme_classic() +
  theme(aspect.ratio=1) + theme(legend.position = "none") +
  scale_color_manual(values = myPal) +
  scale_shape_manual(values = myPalShape) +
  labs(color="Intelligibility Type",
       shape = "Intelligibility Type",
       linetype = "Intelligibility Type")

# Creating OT Scatterplot Figure

scatter <- VSA  + disp + autoVSA + Hull + vsd25 + vsd75 +
  patchwork::plot_layout(guides = 'collect',
                         ncol = 3) & theme(legend.position = "bottom")
scatter 

ggsave("Plots/ModelFigure.png", scatter,
       height = 4.5,
       width = 6,
       units = "in",
       scale = 1.1)

Filtering at Different Levels

    Hull_b <- cHull(Formants_PRAAT$F1_b, Formants_PRAAT$F2_b)

Error in Polygon(chull.coords, hole = F) : 
  could not find function "Polygon"

Listener Demographic Information

ListenerDemo <- Listeners %>%
  furniture::table1(age, gender, race, ethnicity)

ListenerDemo

Speaker Demographics

SpeakerDemo <- AcousticData %>%
  dplyr::select(c(Speaker, Sex, Etiology))

Ages <- rio::import("Prepped Data/Speaker Ages.xlsx")

SpeakerDemo <- full_join(SpeakerDemo, Ages, by = "Speaker")

SpeakerDemoInfo <- SpeakerDemo %>%
  furniture::table1(Sex, Etiology, Age, na.rm = F)

SpeakerDemoInfo

SpeakerDemo %>%
  dplyr::summarize(mean_age = mean(Age, na.rm = T), age_sd = sd(Age, na.rm = T), age_range = range(Age, na.rm = T))